function r = laborassignment(sbar,stilde,xlow,xup,p)

% check for consistent automation thresholds and correct if necessary
if xlow > xup
    xup = xlow; % set xup to xlow if the supplied xup is below xlow
    warning('The supplied xup is smaller than xlow. Continue with xup reset to xlow.');
end

% compute bottom branches of log chi (defined as lambda*log(w)-log(y)) and
% x
if stilde>sbar
    
    % find chi_stilde^- (defined as limit of w_s^lambda/y as s->stilde from below)
    chistildeminus = fzero(@(chistart) terminalvalue(chistart,[stilde stilde/2+sbar/2 sbar],xlow,p),[xlow*p.stol,10/p.stol],p.optionsfzero);
    
    % use chi_stilde^- to obtain log(chi) 
    % and assignment function
    [sgridbottom, logchixbottom] = ode45(@(s,logchix) diffeqchi(s,logchix,p),[stilde sbar],[log(chistildeminus) xlow]);
    
    % extract log chi and assignment on [0,stilde]
    sbottom = flip(sgridbottom);
    logchibottom = flip(logchixbottom(1:end,1));
    xbottom = flip(logchixbottom(1:end,2));
    
    % compute labor supply on grid (normalized to one in baseline, so only needed
    % when doing comparative statics with respect to labor supply)
    if p.labor == "new"
        laborbottom = labornew(sbottom);
    end
    
end

% compute top branches of log chi and x
if stilde<1
    
    %find chi_stilde^+ (defined as limit of w_s^lambda/y as s->stilde from above)
    chistildeplus = fzero(@(chistart) terminalvalue(chistart,[stilde stilde/2+1/2 1],xup,p)-1,[(1-xup)*p.stol,10/p.stol],p.optionsfzero);
    
    % use chi_stilde^- to obtain log(chi) (defined as lambda*log(w)-log(y))
    % and assignment function
    [sgridtop, logchixtop] = ode45(@(s,logchix) diffeqchi(s,logchix,p),[stilde 1],[log(chistildeplus) xup]);
    
    % extract log chi and assignment on [0,stilde]
    stop = sgridtop;
    logchitop = logchixtop(1:end,1);
    xtop = logchixtop(1:end,2);
    
    % compute labor supply on grid (normalized to one in baseline, so only needed when doing comparative statics with respect to labor supply)
    if p.labor == "new"
        labortop = labornew(stop);
    end
    
end

% compute output, wages, and collect results
if stilde>sbar && stilde<1
    
    % compute gamma_k
    gammak = integral(@(x) exp((p.lambda-1)*(logcapprod(x,p)-p.q)), xlow, xup);
    
    % check for positive, finite output
    if p.lambda<1 && (gammak*exp((p.lambda-1)*p.q))>=1 % zero output
        
        r.sgrid = [sbottom; stop];
        r.sbottom = sbottom;
        r.stop = stop;
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = zeros(length(r.sgrid),1);
        r.wbottom = zeros(length(sbottom),1);
        r.wtop = zeros(length(stop),1);
        r.alphak = [];
        r.y = 0;
        r.c = 0;
        
        return
        
    elseif p.lambda>1 && (gammak*exp((p.lambda-1)*p.q))>=1 % infinite output
        
        r.sgrid = [sbottom; stop];
        r.sbottom = sbottom;
        r.stop = stop;
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = inf*ones(length(r.sgrid),1);
        r.wbottom = inf*ones(length(sbottom),1);
        r.wtop = inf*ones(length(stop),1);
        r.alphak = [];
        r.y = inf;
        r.c = inf;
        
        return
        
    end
    
    % compute output
    if p.labor == "new" % using new labor supply vector as given by function labornew.m (for comp statics wrt labor supply)
        y = ((trapz(sbottom,exp(logchibottom./p.lambda).*laborbottom)+trapz(stop,exp(logchitop./p.lambda).*labortop))/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    else % using baseline labor supply, normalized to one
        y = ((trapz(sbottom,exp(logchibottom./p.lambda))+trapz(stop,exp(logchitop./p.lambda)))/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    end

    
    % compute wages
    wbottom = exp(logchibottom).^(1/p.lambda)*y^(1/p.lambda);
    wtop = exp(logchitop).^(1/p.lambda)*y^(1/p.lambda);
    
    % compute capital share
    alphak = gammak*exp(p.q*(p.lambda-1));
    
    % merge bottom and top branches of sgrid, w, x
    sgrid = [sbottom; stop];
    w = [wbottom; wtop];
    x = [xbottom; xtop];
    
    % collect results
    r.sgrid = sgrid;
    r.sbottom = sbottom;
    r.stop = stop;
    r.x = x;
    r.xbottom = xbottom;
    r.xtop = xtop;
    r.w = w;
    r.wbottom = wbottom;
    r.wtop = wtop;
    r.alphak = alphak;
    r.y = y;
    r.c = (1-alphak)*y;
    
elseif stilde==sbar
    
    % compute gamma_k
    gammak = integral(@(x) exp((p.lambda-1)*(logcapprod(x,p)-p.q)), 0, xup);
    
    % check for positive, finite output
    if p.lambda<1 && (gammak*exp((p.lambda-1)*p.q))>=1 % zero output
        
        r.sgrid = stop;
        r.sbottom = [];
        r.stop = stop;
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = zeros(length(r.sgrid),1);
        r.wbottom = zeros(length(r.sgrid),1);
        r.wtop = zeros(length(stop),1);
        r.alphak = [];
        r.y = 0;
        r.c = 0;
        
        return
        
    elseif p.lambda>1 && (gammak*exp((p.lambda-1)*p.q))>=1 % infinite output
        
        r.sgrid = stop;
        r.sbottom = [];
        r.stop = stop;
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = inf*ones(length(r.sgrid),1);
        r.wbottom = [];
        r.wtop = inf*ones(length(stop),1);
        r.alphak = [];
        r.y = inf;
        r.c = inf;
        
        return
        
    end
    
    % compute output
    if p.labor == "new" % using new labor supply vector given by labornew.m (for comparative statics wrt labor supply)
        y = (trapz(stop,exp(logchitop./p.lambda).*labortop)/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    else % using baseline labor supply normalized to one
        y = (trapz(stop,exp(logchitop./p.lambda))/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    end
    
    % compute wages
    w = exp(logchitop).^(1/p.lambda)*y^(1/p.lambda);
    
    % compute capital share
    alphak = gammak*exp(p.q*(p.lambda-1));
    
    % collect results
    r.sgrid = stop;
    r.sbottom = [];
    r.stop = stop;
    r.x = xtop;
    r.xbottom = [];
    r.xtop = xtop;
    r.w = w;
    r.wbottom = [];
    r.wtop = w;
    r.alphak = alphak;
    r.y = y;
    r.c = (1-alphak)*y;
    
elseif stilde==1
    
    % compute gamma_k
    gammak = integral(@(x) exp((p.lambda-1)*(logcapprod(x,p)-p.q)), xlow, 1);
    
    % check for positive, finite output
    if p.lambda<1 && (gammak*exp((p.lambda-1)*p.q))>=1 % zero output
        
        r.sgrid = sbottom;
        r.sbottom = sbottom;
        r.stop = [];
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = zeros(length(r.sgrid),1);
        r.wbottom = zeros(length(sbottom),1);
        r.wtop = [];
        r.alphak = [];
        r.y = 0;
        r.c = 0;
        
        return
        
    elseif p.lambda>1 && (gammak*exp((p.lambda-1)*p.q))>=1 % infinite output
        
        r.sgrid = sbottom;
        r.sbottom = sbottom;
        r.stop = [];
        r.x = [];
        r.xbottom = [];
        r.xtop = [];
        r.w = inf*ones(length(r.sgrid),1);
        r.wbottom = inf*ones(length(sbottom),1);
        r.wtop = [];
        r.alphak = [];
        r.y = inf;
        r.c = inf;
        
        return
        
    end
    
    % compute output
    if p.labor == "new" % using new labor supply given by labornew.m, for comp statics wrt labor supply
        y = (trapz(sbottom,exp(logchibottom./p.lambda).*laborbottom)/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    else % using baseline labor supply, normalized to one
        y = (trapz(sbottom,exp(logchibottom./p.lambda))/(1-gammak*exp((p.lambda-1)*p.q)))^(p.lambda/(p.lambda-1));
    end
    
    % compute wages
    w = exp(logchibottom).^(1/p.lambda)*y^(1/p.lambda);
    
    % compute capital share
    alphak = gammak*exp(p.q*(p.lambda-1));
    
    % collect results
    r.sgrid = sbottom;
    r.sbottom = sbottom;
    r.stop = [];
    r.x = xbottom;
    r.xbottom = xbottom;
    r.xtop = [];
    r.w = w;
    r.wbottom = w;
    r.wtop = [];
    r.alphak = alphak;
    r.y = y;
    r.c = (1-alphak)*y;
    
end

end